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We apply lattice Boltzmann methods to study the segregation of binary fluid mixtures under 
oscillatory shear flow in two dimensions. The algorithm allows to simulate systems whose dynamics 
is described by the Navier-Stokes and the convection-diffusion equations. The interplay between 
several time scales produces a rich and complex phenomenology. We investigate the effects of 
different oscillation frequencies and viscosities on the morphology of the phase separating domains. 
We find that at high frequencies the evolution is almost isotropic with growth exponents 2/3 and 1/3 
| in the inertial (low viscosity) and diffusive (high viscosity) regimes, respectively. When the period 

of the applied shear flow becomes of the same order of the relaxation time Tr of the shear velocity 
profile, anisotropic effects are clearly observable. In correspondence with non-linear patterns for 
the velocity profiles, we find configurations where lamellar order close to the walls coexists with 
isotropic domains in the middle of the system. For particular values of frequency and viscosity it 
can also happen that the convective effects induced by the oscillations cause an interruption or a 
. slowing of the segregation process, as found in some experiments. Finally, at very low frequencies, 

the morphology of domains is characterized by lamellar order everywhere in the system resembling 
what happens in the case with steady shear. 
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I. INTRODUCTION 



The process of segregation in fluid mixtures is greatly affected by the presence of applied flows . The domains 
■ of the separating phases generally grow with anisotropic patterns that reflect the profile of the velocity field. In the 
case of polymer blends subjected to a simple shear flow string-like structures are observed aligned along the flow 
direction yQ. In more complex systems like diblock copolymer melts with lamellar order 0] the question of the 
orientation of the interfaces is still a debated problem |]J7]|| . The presence of the flow has also less expected and 
obvious consequences. For example, in phase separation of binary fluids, while in absence of flow the size of domains 
is distributed around a single average length-scale |J, when a shear flow is applied, two typical lengths are observed 
for each direction Il0yil|jl^ , |l3|| . Another peculiar case is that of lamellar sheared systems where the symmetry of 
dynamical scaling K generally holding in ordering systems, is foreseen to be violated Related to the presence of 
the flow is also the behavior of the stress response. The non-linear character of this response reflects the dynamical 
I i evolution of the domain pattern and is of fundamental importance in many technological applications . 

In this paper we study through numerical simulations the behavior of a symmetric phase separating binary mixture 
subjected to an oscillatory shear flow jjj. Experiments on this system show that in some cases the growth of the 
domains is interrupted for frequencies smaller than some critical value |l5| , |l6| ] while in other cases domains are observed 
to grow on time scales much longer than the period of a single oscillation jL7). Available simulations of this system 
p3| , ^8|Jl9| did not take into account the role of hydrodynamics or the existence of a finite time required to set a linear 
shear profile in a fluid system. This time, which is inversely proportional to the viscosity in a simple fluid [ po[ , has a 
very important role in the case of oscillatory shear. For example, for sufficiently high frequencies, as it will be shown 
later, this time is longer than the oscillation period and the linear profile will be never set in the system. This will 
have relevant consequences for the evolution and the morphology of the domains of the separating phases. Actually, 
a systematic study of the dependence of the growth properties on the frequency of oscillations has not been done 
in previous simulations also for the simplest cases without hydrodynamics. One has also to observe that, differently 
from the case with steady shear, due to the fact that the average strain is zero in one period of shearing, the effects 
of oscillatory shear on the morphology of domains are not always easy to understand intuitively. 

The effects discussed above can be properly described only considering the full hydrodynamic equations for binary 
mixtures. We have used Lattice Boltzmann Method (LBM) |2l],^,^3j to simulate the convection-diffusion and Navier- 
Stokes equations for a binary fluid. We have introduced in the lattice Boltzmann scheme appropriate boundary 
conditions for a shear flow and we have run our simulations systematically changing the frequency of oscillations for 
a limited set of values of the viscosity. We have considered the two-dimensional case which is also useful for the 
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comprehension of the three-dimensional case and has the advantage of being less demanding from the computational 
point of view. 

Lattice Boltzmann simulations are based on discretized Boltzmann equations for a set of distribution functions 
related to the fluid densities and velocity. The densities and momenta are conserved at each step of the simulation 
thus approximating the behavior of the hydrodynamic equations for the fluid. Lattice Boltzmann methods have been 
found to be very convenient for simulating quasi-incompressible fluids on very long time scales, as it is needed in 
phase separation problems 24 2l|. Another advantage of the LBM is that, in the implementation of the method used 
in this paper |p6f , a free energy can be introduced such that the fluid relaxes to the equilibrium state determined by 
this free energy. This allows to know accurately the equilibrium properties of the coexisting phases whose growth 
dynamics is under study. 

Before presenting our results it is useful to summarize the known behavior of two-dimensional binary mixtures 
quenched into the ordered homogeneous phases in absence of flow. Once domains of the two phases are well established 
experimental and numerical data show that the growth is self-similar with the typical size of domains scaling with time 
as R(t) ~ t a |J. The growth exponent a depends on the physical mechanism responsible for the phase separation. 
A simple scaling analysis of the Navier-Stokes and of the convection-diffusion equations shows that three regimes 
can be found corresponding to the role played by hydrodynamic degrees of freedom |27],2^]. At high viscosity the 
domain growth is governed by a diffusive mechanism and the growth exponent is a — 1/3 29| . When hydrodynamics 
becomes relevant, the laws Rjt) ~ t or R(t) ~ i 2 / 3 are expected depending on whether viscous forces or inertial effects 
dominate, respectively |27|]2q |. In real systems, however, the situation is more complex. The physical mechanism 
responsible for the viscous growth is not operating in the two-dimensional case p0[ | and, indeed, this regime has never 
been observed in simulations 0] . 

The effects of a steady shear flow on the growth laws previously discussed have been considered in many papers. 
In the diffusive case, analytical calculations based on a self-consistent approximation show that the typical size of 
domains should grow in the direction normal to the flow as in the case without shear while the growth exponent 
in the flow direction is equal to that in the transverse direction augmented of one |I~i| , |3"l| ] . This result cannot be 
easily checked by numerical simulations due to the presence of finite size effects that become very soon relevant in 
the direction of the flow affecting the value of the exponents J^L The full problem with the Navier-Stokes and the 
convection-diffusion equations has been considered in Refs. p5[ 32| , [33" 34 3^, 3(| but also in this case reliable results 
for the growth exponents are not yet available. Actually, the true existence of an asymptotic growth regime in 
alternative with a stationary state with a finite transverse size is a question still to be clarified |]37| , p5f . On the other 
hand, morphological properties are reasonably well understood. Domains are stretched by the flow and this induces 
a coagulation of domains in the flow direction but also ruptures in the bicontinuous network |l^ , |38|| . As a result, 
domains assume the typical string-like character with the already mentioned complication that the size of domains is 
distributed around two typical length-scales for each direction. 

The lattice Boltzmann scheme used in this paper is described in the following Section. Due to small variations 
with respect to previous LBM |35j, details on the method and on the implementation are given for convenience of 
the reader. The rest of the paper is divided as follows. In Section III we illustrate our strategy for the choice of the 
parameters used in the simulations; we also discuss the relevant time scales for the system considered. In Section IV 
we show results of simulations where the growth is dominated by inertial effects while the diffusive case is treated in 
Section V. In Section VI we consider the behavior of the shear stress and we draw our conclusions in Section VII. 



II. THE MODEL 

Our simulations are based on the lattice Boltzmann scheme developed by Orlandini et al. |3{J and Swift et al. jlO| . 
In this scheme the equilibrium properties of the system can be controlled by introducing a free energy which enters 
properly into the lattice Boltzmann model. 



A. The lattice Boltzmann scheme 



The lattice Boltzmann scheme is based on the D2Q9 lattice: A square lattice is used in which each site is connected 
with nearest and next-to-nearest neighbors. The horizontal and vertical links have length Ax, the diagonal links 
\/2Ax, Ax being the space step. Two sets of distribution functions fi(r,t) and gi(r,t) are defined on each lattice site 
r at each time t. Each of them is associated with a velocity vector e^. Defined At as the simulation time step, the 
quantities e^At are constrained to be lattice vectors so that |ej| = Ax/ At = c for i=l (East direction), 2 (North), 
3 (West), 4 (South) and |e;| = \/2c for i=5 (North East), 6 (North West), 7 (South West), 8 (South East). Two 
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functions /o(r,i) and go(r,t), corresponding to the distribution components that do not propagate (eo = 0), are also 
taken into account. They evolve during the time step At according to a single relaxation-time Boltzmann equation 
©[§: 

Mr + ei At,t + At) - Mr,t) = --[/i(r,i) - f?(r,t)], (1) 

T 

9i (r + e,At, t + At)- 9l {r, t) = - — [ gi (r, t) - gf (r, t)\, (2) 

where r and r v are independent relaxation parameters, f^ q (r,t) and g^ q (r,t) are local equilibrium distribution func- 
tions. The distribution functions are related to the total density n, to the fluid momentum nv and to the density 
difference ip through 

n = ^2fi, nv = fi^h <f = 9i ' ^ 

i i i 

These quantities are locally conserved in any collision process and, therefore, we require that the local equilibrium 
distribution functions fulfill the equations 

i i 

YM q -9i) = ^^T.^ q = ^ ( 4 ) 

i i 

E(/r-/i)e, = 0^^/fe i= nv 

i i 

Following Refs. |3^,^], the higher moments of the local equilibrium distribution functions are defined so that the con- 
tinuum equations pertinent to a binary fluid mixture can be obtained and the equilibrium thermodynamic properties 
of the system can be controlled. We define 

^ /^ 9 e iQ e 4/ 3 = c 2 P af 3 + nv a V0 , (5) 

i 

Y 9t q ^ a = <pv a , (6) 
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'eiaeip = c 2 TA^5 a( 3 + tpv a vp . (7) 



where P a p is the pressure tensor, A/x is the chemical potential difference between the two fluids and T is a coefficient 
related to the mobility of the fluid. We want to stress that we are considering a mixture with the two fluids having 
the same mechanical properties and, in particular, the same viscosity. The constraint (|^) expresses the fact that the 
two fluids have the same velocity. The local equilibrium distribution functions can be expressed as an expansion at 
the second order in the velocity v 

f^ = A + C v 2 

jl q = Ai + Biv a e ia + Civ 2 + D I v a vpe ia e i p + G^ a pe ia eip i = l,2,3,4 (8) 
f- q = An + B n v a e ia + C n v 2 + D H v a vpe ia eip + G u ,<xpZia.e-ip i = 5, 6, 7, 8 

and similarly for the g eq ', i = 0, 8. The relations fl)-(|5|) can be used to fix the coefficients of these expansions. A 
suitable choice of the coefficients in the expansions (J8p is 

A = n-20A II , A I = 4A II , A n = (9) 
Bj = 4B Hl B u = ^2 ( 10 ) 
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C = C l = AC U, Cjl = (11) 

Di = 4D n , D n = £j (12) 

n —ah n _ Pgf3 ~ \P<j<jQ~a(3 , . 

^1,0/3—^11,0/3, ^II, a p — ' 

The expansion coefficients for the g^ q can be got from the previous ones with the formal substitutions n — > tp and 

Pa/3 — * rA/i<5 Q ^. The quantities P a /3 and A/i, which appear in the coefficients of the equilibrium distribution functions, 
can be calculated from a suitable free energy. 

B. The equilibrium properties 



The free-energy functional used in the present study is 



T = \ dr 



-nlnn + -ip 2 + -ip 4 + -(V(p) 2 
3 2 y V 2 V yj 



(14) 



The term in n gives rise to a positive background pressure and does not affect the phase behavior; it is required in 
the lattice Boltzmann approach, as it will be seen later. The terms in <p in the free-energy density f(n, ip) correspond 
to the usual Ginzburg-Landau free energy typically used in studies of phase separation |28| . The polynomial terms 
are related to the bulk properties of the fluid. While the parameter b is always positive, the sign of a distinguishes 
between a disordered (a > 0) and a segregated mixture (a < 0) where two pure phases with ip = ±y/—a/b coexist. Wc 
will consider quenches into the coexistence region with a < and b — — a so that the equilibrium values for the order 
parameter are <p = ±1. The gradient term is related to the interfacial properties. The equilibrium profile between the 

two coexisting bulk phases is (fi(x) = tanh \ — x giving a surface tension 

V 2/t 



-V^2aH (15) 



and an interfacial width 



e = 2j — ■ (16) 
V —a 

The thermodynamic properties of the fluid follow from the free energy (|l4|). The chemical potential difference 
between the two fluids is given by 

A^i = — = aip + bip 3 - kVV (17) 
5(p 

The pressure is a tensor P a p since interfaces in the fluid can exert non- isotropic forces |44|] . The diagonal part p can 
be calculated from ( |l4| ) by using thermodynamics relations: 

ST 5T 

Po = n — +ip }{n, ip) 

on dip 

= \n+ ^ 2 + ^ - ^ (W) ~ \ (V^) 2 (18) 

In deriving the pressure tensor P a p, one has to ensure that P a p obeys the condition of mechanical equilibrium 
d a P a (3 — JH]. A suitable choice is 

P a /3 = Po$ a p + nd a ipdpip. (19) 

The presence of the term depending on n in the free-energy density allows to recover the known results of the D2Q9 
lattice Boltzmann model for a single fluid. Indeed, when a = b = n = ip = 0, the expansion coefficients (p|)-(|lq) 
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become those of the D2Q9 model and the pressure tensor (|l9|) reduces to the scalar p = (c 2 /3)n, where we have 
also included the factor c 2 appearing in the r.h.s. of Eq. (||). This is the pressure for an ideal gas with speed of sound 
c s = cj\f% p6| . Let us observe that the value of the numerical factor in front of nhin in the free energy depends on 
the topology of the lattice and the spatial dimensions [fl7) . 

It has been shown in Refs. J39pc| ], using a Chapman-Enskog expansion jjfifl , that the above described lattice 
Boltzmann scheme simulates at second order in At the continuity, the quasi-incompressible Navier-Stokes and the 
convection-diffusion equations with the kinematic viscosity v and the macroscopic mobility O given by 

v = Atj(r-±), 9 = rAte 2 (r v -i). (20) 



C. The shear boundary condition 



In order to enforce a shear flow on the system, we have used the following scheme. We assume that the shear 
flow is directed along the horizontal direction. Boundary walls are placed on the upper and lower rows of sites. Let 
us consider the upper wall (similar considerations apply to the lower wall). After the propagation the distribution 
functions fo(t), fi(t), hit), hit), hit) and hit) are known on each site, while hit), hit), f$(t) are not. One uses 

Eqs. (||) to determine them as well as n. Requiring that the wall velocities w Xt t = 70 — ^ — cos(2ir ft), w y .t — are 
imposed to the fluid, we can write 

fr(t) + h(t) + f a {t) =n- [/ (f) + h{t) + f^t) + hit) + /„(*) + hit)] 

hit) - f 7 {t) = n 7o ^ cos(2^/t) - [Mt) - Mt) + hit) - f^t)} (21) 

frit) + hit) + hit) = hit) + hit) + h{t) 

where L is the lattice size, 70 is the amplitude of the shear rate and / is the frequency of the oscillatory shear. 
Consistency of Eqs. (|2l]) gives 

n = hit) + hit) + hit) + 2 [hit) + hit) + hit)} (22) 

The system of Eqs. (pi|) reduces to two equations with three unknown variables. To close the system of equations 
the bounce-back rule [49|j50[ | is adopted for the distribution functions normal to the boundary. This corresponds to 
impose that f^it) = hit)- In order to preserve correctly mass conservation we add a further constraint. Mass will be 
conserved if the total density n on each site is equal to the quantity n given by the sum 

hit,t-At) = hit ~ At) + hit - At) + f 2 {t - At) + hit - At) 

+/i(t) + fs(t) + hit) + hit) + hit). (23) 

where quantities at time (i — At) refer to the previous time step and have not been propagated over the lattice. In 
order to impose the constraint that on all the boundary sites n = n, we have to introduce an independent variable 
in the system of equations. We have chosen hit) since it does not propagate |34). The solutions of the system of 
Eqs. ( pl| ) and n = h are 

/o(r, t) = n - [h (r, t) + / 3 (r, t)] - 2 [/ 2 (r, t) + h (r, t) + f 6 (r , t)} 
hir,t) = h(r,t) 

h (r , t) = /„ (r , t) - X - [f x (r, f ) - / 3 (r, t)\ + \ n 7o ^ cos(2^/<) (24) 
/ 7 (r,i) = / 5 (r,i) + 1 [A(r,<) - / 3 (r,t)] - 1 n 7o cos(2^/<) 

A similar treatment is required for the giir,t). With this choice the proper momentum at the boundary is achieved. 
At this point the collision step is applied to all sites, including the boundary ones. By this procedure, once the system 
has been initialized, the application of the propagation and collision steps goes on preserving mass and momentum 
conservation and implementing the correct velocity values on the boundaries, as it has also been verified numerically 
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Finally, we also require that the two fluids, which, as already stated, are assumed to have the same dynamic and 
static properties, have a neutral wetting with walls. This can be enforced at each time step by the condition 



m-V(y5 = (25) 

where m is a unit vector normal to the wall |5l|,p2[ . This corresponds to fix the gradient of the density (p on the 
walls so that the angle formed by the interfaces between the two fluids with the walls stays at a constant value of w/2 
radians. This completes the description of the model used in the present work. 

III. PARAMETER SELECTION AND RELEVANT TIME SCALES 

We have studied the effects of the applied flow by changing the frequency of the oscillations at fixed values of the 
parameters a, b, k, t, V. This has been systematically done for the cases reported in Table |. We have kept the ratio 
n/a fixed in such a way that the interfacial width is always of about 3 lattice spacings. We fixed t v — 1 with T 
controlling the value of the macroscopic mobility (^Oj). The amplitude of the shear rate is equal in all runs to the 
value 70 = 0.005. The size of the lattice, if not otherwise stated, is L = 256. 

In absence of shear, as observed in |^5|,[53| , the behavior at late times of a viscous phase separating binary mixture can 
be described in terms of adimensional temporal and spatial quantities. Indeed, from the set of macroscopic parameters 
n,a,iy it is possible to define only one unit of length (To = nv 2 ja) and one of time (To = n 2 v /a 2 ) II Then, m a 
regime of dynamical scaling with the size of domains distributed around a single typical length, when diffusion is 
negligible, it is possible to build up only one spatial and one temporal adimensional variables. These variables have 



been used in comparing results of simulations performed with different parameters and methods [|25 53 1. 

The same use of adimensional variables can be done in the case of very high viscosity when the evolution equations 
reduce to a convection-diffusion equation. To = £ 3 /(6cr) IS the time scale for diffusion f2g| . 

In the case with oscillatory shear the situation is more complex because various temporal scales can be defined 
and in general one does not expect dynamical scaling. Due to these reasons we have preferred to show our results 
in terms of original quantities. However, it remains useful to consider the relevant time scales because they can give 
informations on the physical mechanisms responsible of phase separation and on the role of the applied flow. 

We have already defined the quantities To and To which do not depend on the applied flow. Then there are the 
time Tg — 7q ; related to the amplitude of the oscillating shear, and the period T = f^ 1 of a single oscillation. 
The ratio Tq/Ts can be interpreted as an indicator of the relative relevance of hydrodynamic and imposed velocities. 
Finally we consider the quantity Tr = L 2 /{vtt 2 ) which is the leading contribution to the relaxation time for a steady 
shear profile in a simple fluid [Q . This time can be taken as indicative of the relaxational velocity processes also in 
the case of a phase separating binary mixture [ |34[ . In Table [j] the time scales corresponding to the sets of parameters 
indicated are also reported. 

The relevance of the ratio Tr/T for the problem considered in this paper appears clearly from Fig. 1. Here the 
horizontal velocity v x profiles for the vertical cross section in the middle of the system are reported at a late stage 
of simulations for the cases 1 and 3 of Table | and for two different frequencies. We checked that these results are 
independent of the particular vertical line considered. For each case the velocity profile is plotted at the quarters of 
one period. For the set of parameters 1 with / = 10 -3 the horizontal velocity induced by the shear is very small 
in the bulk of the system and comparable with the average of the modulus of the vertical velocity v y . With the 
same parameters, when the frequency decreases to / = 10~ 4 , the complex non-monotonic behavior of the horizontal 
velocity is more evident. v x is much larger than hydrodynamic velocities in a relevant portion of the system close to 
the walls where it also assumes opposite directions to those imposed by the walls. This peculiar pattern of the velocity 
profiles has consequences for the behavior of the stress, as it will be seen later. When the viscosity becomes higher, 
like for the other set of parameters used in Fig. 1, the relaxation time Tr decreases. Then, for the same frequencies, 
the velocity in the bulk of the system is larger and almost triangular oscillating profiles can be observed at / = 10~ 4 
. We have also checked that each set of four profiles of Fig. 1 is typical for the parameters considered in the sense 
that only small variations in the pattern of these profiles, probably due to the evolving interface configuration, can 
be observed during the simulation. Figure 2 shows the velocity profile at the same phase of four different periods at 
initial and late stages of the simulation with set 1 of parameters at / = 10 -3 . Small quantitative changes can be 
observed in the region close to the walls while the general shape of the profile always remains the same. 



1 These expressions do not depend on the spatial dimension. 
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Finally we consider the question of the stability of the lattice Boltzmann scheme used in this work. As observed in 
|f33f , this lattice Boltzmann scheme is intrinsically unstable. At unpredictable times of the simulation pressure waves 
grow up indefinitely in very few iterations making not possible the continuation of the simulation. As expected, we 
saw that this problem becomes more serious when r tends to the limit 1/2. The problem of stability is connected to 
the one of guaranteeing as much as possible the incompressibility of the fluid. Compressibility errors, which go like 
(v/c s ) 2 [p4| , can be reduced by either increasing c s , which would require to reduce At, or decreasing the magnitude of 
a, 6, k [|53|]. We have followed a mixed strategy by keeping Ax = 1 and changing the values of At as reported in Table |. 
In this way we were able to run simulations long enough to study the phase separation of binary mixtures in different 
growth regimes implicitly assuming that the evolution of the system is not affected by the possible occurrence of the 
numerical instability. A comment to the results of Ref. |53| is that the introduction of walls for the shear boundary 
conditions makes worst the stability properties of the LB scheme. 



IV. INERTIAL ORDERING 



In two-dimensional quiescent systems, as discussed in the introduction, tweugrowth regimes with different power 
law behaviors for the average size of domains R(t) have been clearly identified □. In this Section we will consider the 
effects of the oscillatory shear on the case of phase separation driven by inertial growth. We will mostly refer to the 
case 1 of Table | for which, in absence of flow, the behavior of R(t) is shown in Fig. 3. The quantity R(t) has been 
calculated as the first momentum of the structure factor, that is 

fdk C(k,t) , . 

*® = J dk k c{k!t) (26) 

where C(k,t) is the spherical average of the structure factor 

= (¥>(MM-M)> (27) 

and (■) is the average over different histories. We found a growth exponent a = 0.62 i; the small discrepancy from 
the expected value a = 2/3 typical for the inertial growth can be attributed to a small violation of dynamical scaling 

n. 

The effects of the oscillatory shear on the growth of the domain size for the case 1 of Table | can be seen in the 
panel of Fig. 4 which summarizes our results for a range of frequencies from / = 10~ 3 up to / = 5 • 1CP 6 . We measure 
the spherical average R{t), the domain size in the flow direction 

Jfc(k,t) 

J dk\k x \C(k,t) 

and the analogue R y for the shear direction. The value / = 10 -3 is the highest frequency where an anisotropic 
behavior can be observed. R and R x evolve with an exponent which is equal to 2/3; the change of the slope of log 10 R 
and log 10 R x at log 10 t ~ 4.2 is due to finite size effects which are more relevant in the direction of the flow. The 
behavior of R y departs from the power law i 2 / 3 at log 10 t ~ 3.5. This anisotropy is partially due to the presence of the 
walls which, even without flow, can inhibit the growth of domains in the vertical dimension at the bottom and at the 
top of the system. The morphology of the domains is also influenced by the shear velocity which is larger compared 
with the hydrodynamic velocities in the region close to the walls - see Fig. 1. The slightly anisotropic evolution at 
/ = 10 -3 can be illustrated from the two configurations shown in Fig. 5 with the corresponding structure factors. 
The circular shape of the structure factor at t = 1650, whose radius is of the order of the inverse of R(t), reflects 
the isotropic configuration of the concentration field at this time. At t = 5900 the slight prevalence of the peaks at 
k x = corresponds to the presence of a recognizable amount of domains aligned with the flow close to the walls. 
In the following of this simulation the almost isotropic character is conserved as it has been checked looking at the 
configurations. 



2 Growth regimes with exponent 1/2 have also been reported; their existence at asymptotic times is still under debate. For a 
discussion see Ref. p4| . 
3 For a more accurate measure of this exponent we have used larger lattices with L = 512. 
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At the frequency / = 10~ 4 the morphology of domains is more affected by the applied flow due to a bigger region 
of the system where the horizontal velocity is significantly larger than the typical hydrodynamic velocities measured 
along the vertical direction - see Fig. 1. As a consequence lamellae can be observed close to the walls while the growth 
keeps a more isotropic character in the middle of the system. Figure 6, at the frequency / = 2 • gives an example 
of this behavior with the evolution of the system shown for a whole period. Larger and more spherical domains can 
be observed in the middle of the system while, close to the walls, thin domains follow the direction of the flow and 
are subjected to a larger number of recombination and breakup processes. In this case and also at / = 5 • 10~ 6 it is 
not possible to speak of dynamical scaling since domains are distributed on different scales. However, it is worth to 
observe that the quantities R and R y follow for a large interval of the evolution a power law behavior with exponent 
2/3. 

When the viscosity becomes larger, as for example in the cases 2,4,5 of Table |. the evolution of the system without 
flow still corresponds to the inertial regime but the resulting shorter relaxation time Tr makes the presence of the 
flow in the bulk more effective with relevant consequences for the kinetics of phase separation. In particular we 
have observed that the oscillatory flow can cause the interruption of the segregation process. This phenomenon, also 
reported in experiments ll^Jl^l , has been found in our simulations at different viscosities and in different growth 
regimes. An example of this flow- induced interruption of growth is shown in Fig. 7. The growth of R, R x , R y becomes 
very slow at log 10 t ~ 4.8. In Fig. 7 we also show a set of 4 configurations in a period at this time. We observe 
that the terminal regions of domains close to the walls follow with their movement the oscillation of the flow. The 
convection-induced movements inhibit the domain growth due to other mechanisms (diffusion or inertial) and the 
system appears for a certain interval of time in a sort of elastic stationary state. The size of domains where this 
phenomenon is first observed during the phase separation is found to be of the order of the average deformation 

/ dt dy~/oy coa(2n ft) — 70 T in all cases considered. 

L - 1 Jo Jo 8?r 

V. DIFFUSIVE GROWTH 

In this Section we consider the case where diffusion is the physical mechanism mainly responsible for phase separa- 
tion. We will consider the set 3 of parameters of Table E the corresponding behavior of R(t) in quiescent conditions 
is shown in Fig. 3 with the value of the growth exponent given by a — 0.35. 

As in the previous cases also here the growth becomes more anisotropic when the frequency / decreases. However, 
due to a higher value of the viscosity, Tr is smaller and the effects of shear convection are more pronounced and 
observable already for / = 1CP 2 . Indeed, at this frequency, and also at / = 1CP 3 , as it can be seen in Fig. 8, at late 
times in the simulation, R x and R grow faster than i 1 / 3 with an effective exponent a = 0.39 at / = 1 0~ 2 and a = 0.54 
at / = 10~ 3 for R(t). This behavior can be understood by looking at the configurations of the concentration field. In 
Fig. 9 it is shown an example at the time t = 2250 for / = 10~ 3 . Two different phases can be seen to coexist: lamellar 
ordered domains aligned with the flow close to the walls and the usual isotropic pattern of phase separation in the 
middle of the system. This coexistence is reflected in the shape of the structure factor which is circular with two peaks 
at k x = corresponding to the horizontal lamellar domains. Then, as in diffusive phase separation with steady shear, 
striped domains almost aligned with the flow grow in the flow direction faster than in the other directions with an 
exponent larger than 1/3. In the case of Fig. 9, the effective exponents for R x and R, which are quantities averaged 
over the whole system, will depend on the ratio between the volumes of the two coexisting phases. 

By decreasing the frequency, the difference in the behavior of R x and R y becomes more pronounced, as it can be 
seen in Fig. 8 at / = 10 -4 and / = 10~ 5 . The four configurations shown in Fig. 10 for the first period of the evolution 
at / = 10 -4 exhibit elongated domains in the direction of the flow similar to those observed in the case of steady 
shear. This explains the big difference in the values measured for R x and R y in Fig. 8. Of course, also in this case the 
morphology of the domains is strictly related to the behavior of the horizontal velocity profiles shown in Fig. 1. We 
see that an almost regular triangular velocity profile occurs when the ratio between Tr and T is of order 1. Finally, a 
quantitative evaluation of the growth at the late stages of the simulation at / = 10 -4 can be better done by averaging 
the behavior of R x and R y over each period. This gives the exponents a x = 0.65 and a y — 0.19. These results should 
be compared with the case of steady shear where a very slow growth is measured at late times in the shear direction. 

VI. STRESS BEHAVIOR 

In this Section we consider the behavior of the shear stress associated with the deformation of the domain pattern 
induced by the flow. The stress response a xy is calculated as the second momentum of the structure factor: 
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We first discuss a peculiar behavior that we find for the phase of the shear stress. We show in Fig. 11 the time 
evolution of a xy for the same frequencies and parameters of Fig. 4. For convenience the velocity on the upper wall is 
also plotted in Fig. 11. We observe at the frequencies / = 10~ 3 and / = 10~ 4 a phase opposition between the stress 
and the velocity field imposed on the walls of the system. This unusual phase behavior can be explained by looking 
at the velocity profiles of Fig. 1 where we see that the velocity not only never relaxes to the triangular profile but 
also assumes the "wrong" sign in proximity of the walls before the jump to the values imposed on the boundaries. 
Therefore the stress follows the sign of the real velocity field in the system and this explains the "strange" phase 
behavior of the stress. 

The above analysis is confirmed when we look in Fig. 12 at the case / = 10~ 4 corresponding to the set 3 of 
parameters of Table |. In this case almost triangular profiles are obtained - see Fig. 1, and indeed the stress is almost 
in phase with the velocity field on the walls □. 

A more general feature of the behavior of the shear stress can be observed in all the cases shown. We see that the 
initial evolution of the stress is always characterized by the presence of a peak which can be eventually followed by 
other large oscillations. This is clear from the inset of Fig. 11 for / = 10 -3 where two overshoots modulated by small 
oscillations due to the velocity field can be observed. The phenomenon is enhanced at / = 1CP 4 and / = 2 • 10~ 5 
where the time between successive peaks is of the same order of the period of the applied flow and is also present at 
/ = 5 • 10~ 6 confirming that its origin is independent of the oscillations of the flow. Indeed, for the latter mentioned 
frequency, we observe overshoots of a xy while the phase of the applied flow has not changed sign. Similar phenomena 
are also observed in Fig. 12. 

Overshoots of the shear stress have been reported in experiments of phase separation with steady shear [^6| and 
have also been found in simulations |12 1^J3^ |. The phenomenon is interpreted as due to an initial stretching of the 
domains in the direction of the flow to which it corresponds an increase of a xy . At a certain point the deformation 
cannot be sustained by the surface tension and the domains start to break evolving in less stretched configurations. 
This occurs in correspondence of a maximum of a xy . Then the system becomes more isotropic but, after that a 
minimum in <j xy is reached, elastic energy is again stored due to the deformations and another overshoot of a xy can 
be observed. 

The above considerations can be extended also to our case with further complications due to the flow oscillations. 
In particular, at frequencies of the order of the inverse of time between two overshoots of a xy , the relaxation or the 
stretching phenomena discussed above in the case of steady shear are greatly influenced by the flow oscillations. For 
example, in the case / = 2 • 10~ 5 of Fig. 11, the relaxation after the first maximum of a xy occurs in correspondence 
of the change of sign of the imposed velocity and, therefore, one can think that the decreasing of the stress is mainly 
due to the reversed sign of the deformations than to breaking processes in the domain pattern. This results in two 
very well shaped overshoots than those occurring in the case / = 5 ■ 10 -6 which resembles what would occur with 
steady shear. 

VII. CONCLUSIONS 

In this work we have studied the behavior of phase separating binary mixtures subjected to oscillatory shear for 
different viscosities and frequencies of the applied flow. The existence of different physical mechanisms operating in 
phase separation, the anisotropic effects induced by the flow, and a finite relaxation time Tr for the triangular velocity 
profile contribute to giving rise to a very rich phenomenology. In particular the role of viscosity is fundamental because 
both the occurrence of inertial or diffusive growth and the time Tr depend on the viscosity. 

In this complex framework we found that the ratio Tr/T can be used as a convenient parameter for measuring 
the effects of the applied flow. At low viscosity and high frequency, for example, when the ratio Tr/T is larger than 
1, the effects of the shear are limited to a region close to the walls of the system, while in the bulk the growth is 
isotropic as in the case without applied flow. Actually, the most interesting phenomena can be observed when the 
above ratio becomes of order 1. Different phases, with lamellar order in regions close to the walls and isotropically 
oriented domains in the central part of the system, have been observed to coexist and evolve together in this case. 



4 A complete analysis of the viscoelastic properties of phase separating binary mixtures, with the evaluation of the elastic and 
viscous parts of the stress response, is beyond the purposes of this work. 
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For particular values of the frequency and viscosity we have also observed the interruption or slowing of the process 
of segregation: interfaces are convected successively in opposite directions with the net effect of inhibiting any other 
growth mechanism, at least for a significantly large time interval in simulations. Finally, for values of the ratio Tr/T 
much less than 1, the domains grow with lamellar morphology everywhere in the system as in the case of steady flow. 

The question of the existence of power-law behavior for the domain size R x , R y has a definite answer only in some 
cases. For sufficiently high frequencies the power-law behavior of the case without flow, inertial or diffusive depending 
on the viscosity, is generally recovered. At lower frequencies, when different scales for the size of domains are observed 
at the same time in the system, it is not possible to speak of dynamical scaling; however, the quantities R x , R y that 
we measure still give information about the growth behavior. We can note that in the limit of very low frequencies 
we have not found any signal of a stationary state with the size of domains reaching or tending to a finite value. 
Actually we have measured effective growth exponent in the shear direction less than the expected values in absence 
of flow which indicate a continuation of the phase separation also at the late times of our simulations. This could be 
an indication for the case with steady shear suggesting a growth of domains also at asymptotic times. 

We hope that our analysis of the segregation process in binary mixtures under oscillatory shear will stimulate a 
more systematic experimental investigation of these systems. A natural continuation of this study will be its extension 
to the three-dimensional case. 
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TABLE I. Parameters used in simulations and corresponding time scales. The shear time scale Ts is equal in all runs to 200. 



FIG. 1. Horizontal velocity profiles taken at the central vertical line of the system. The profiles are recorded at the beginning, 
half, first and third quarter of the period considered. 

FIG. 2. The evolution of the horizontal profile for the case and at the times reported in the inset. The main frame is a 
magnification of the profiles in the region close to the upper wall. 

FIG. 3. Evolution of the domain size in cases without applied flow with periodic boundary conditions. The size of the lattice 
used in these simulations is L — 512. 

FIG. 4. Evolution of the horizontal, vertical, and spherically averaged size of domains in double logarithmic scale. The 
straight dashed lines have the slope written in the insets. 

FIG. 5. Configuration of domains with the corresponding structure factor at two times in a run with parameter Set 1 
of parameters at the frequency / = 10 . The variables k x (horizontal axis) and k y (vertical axis) vary in the interval 
[— 5/167T, 5/167r] and [— 7r/8, tv/8] at the times t — 1650 and t = 5900, respectively. 



FIG. 6. Four configurations and corresponding velocity profiles in a run with period T = 50000. 

FIG. 7. Four configurations and size domain evolution for a run with parameter Set 4 and / = 10~ 3 . The configurations 
are taken at the four quarters of a period starting at t = 60000. The quantities R m x,Rmy,Rm are the horizontal, vertical, 
spherically averaged domain sizes averaged over each period. 

FIG. 8. Evolution of the horizontal, vertical, and spherically averaged domain size in double logarithmic scale. The straight 
dashed lines have the slope written in the insets. 



FIG. 9. A configuration with corresponding structure factor at t = 2250 in a simulation with parameter Set 3 and / = 10 
In the structure factor plot k x (horizontal axis) and k y (vertical axis) vary in the interval [— 7r/2, 7r/2]. 



FIG. 10. Four configurations in a run with parameter Set 3 and / = 10 



FIG. 11. Behavior of the shear stress compared with the velocity on the upper wall. In the case with frequency / = 10 
the overall behavior is shown in the inset and only the initial evolution is plotted in the main frame. The plot of the velocity 
is translated along the vertical axis for graphical convenience; units are arbitrary for both the quantities. 

FIG. 12. Shear stress and applied flow as in the previous figure, now for cases with parameter Set 3. In the case with 
frequency / = 10~ 3 the overall behavior is shown in the inset while in the main frame the evolution in the time interval 
[10000, 20000] is plotted. The plot of the velocity is translated along the vertical axis for graphical convenience; units are 
arbitrary for both the quantities. 
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